subroutine ReadMol_init( )
use prec
use cell
use general
implicit none
   !----------------------------------------------------------------------------
   integer             :: i, moltype, TCdim = 4, id
   character (len=256) :: molfile
   character (len=1024):: oneline
   logical             :: fexst
   !----------------------------------------------------------------------------
   subnamep = subname
   subname  = 'ReadMol_init'
   !
   if ( allocated( bondsMol  ) ) deallocate( bondsMol  )
   if ( allocated( anglsMol  ) ) deallocate( anglsMol  )
   if ( allocated( dihrsMol  ) ) deallocate( dihrsMol  )
   if ( allocated( attypMol  ) ) deallocate( attypMol  )
   if ( allocated( atposMol  ) ) deallocate( atposMol  )
   if ( allocated( molIDMol  ) ) deallocate( molIDMol  )
   if ( allocated( rotatedMol) ) deallocate( rotatedMol)
   if ( allocated( chargeMol ) ) deallocate( chargeMol )
   if ( allocated( EMass     ) ) deallocate( EMass     )
   if ( allocated( sigma     ) ) deallocate( sigma     )
   if ( allocated( epison    ) ) deallocate( epison    )
   if ( allocated( bondcoef  ) ) deallocate( bondcoef  )
   if ( allocated( anglcoef  ) ) deallocate( anglcoef  )
   if ( allocated( torscoef  ) ) deallocate( torscoef  )
   if ( allocated( DistMin   ) ) deallocate( DistMin   )
   if ( allocated( EName     ) ) deallocate( EName     )
   !
   write(*,'(/,10x,"Please input the file name that carries the molecule info: ", $)')
   ioerr = 1
   do while (ioerr.ne.0)
     read(*, *, iostat=ioerr) molfile
   end do
   inquire(file=molfile, exist=fexst)
   if (.not.fexst) then
      write(*,'(10x,"File:", A, " not found!")') trim(molfile)
      stop
   endif
   open(10, file=molfile, status='old', action='read', iostat=ioerr)
   if (ioerr.ne.0) then
      write(*,'(10x,"Error while opening file:", A, ".")') trim(molfile)
      stop
   endif
   ! First line is the name of the molecule
   read(10, '(A)', iostat=ioerr) molname
   ! 2nd line is the chemical formular of the molecule
   read(10, '(A)', iostat=ioerr) molform
   ! 3rd line: ntypeMol nbndtypMol nangtypMol ndihrsMol
   read(10, *, iostat=ioerr) ntypeMol, nbndtypMol, nangtypMol, ndihtypMol
   if (ioerr.ne.0) then
      write(*,'(10x,"Error while reading # of types info from file:", A, ".")') trim(molfile)
      close(10)
      stop
   endif
   if (ntypeMol.lt.1.or.nbndtypMol.lt.0.or.nangtypMol.lt.0.or.ndihtypMol.lt.0) then
      write(*,'(10x,"Wrong # of types info from file:", A, ": # < 1")') trim(molfile)
      close(10)
      stop
   endif
   ! 4th line: natomMol nbondsMol nanglsMol ndihrsMol
   read(10, *, iostat=ioerr) natomMol, nbondsMol, nanglsMol, ndihrsMol
   if (ioerr.ne.0) then
      write(*,'(10x,"Error while reading # of atoms/bonds/angles/dihedrals info from file:", A, ".")') trim(molfile)
      close(10)
      stop
   endif
   if (natomMol.lt.1.or.nbondsMol.lt.0.or.nanglsMol.lt.0.or.ndihrsMol.lt.0) then
      write(*,'(10x,"Wrong # of atoms/bonds.. info from file:", A, ": # < 1")') trim(molfile)
      close(10)
      stop
   endif
   ! 5th line: density, and optionally the first dimension size of the torsion coefficients, for OPLS, it is 4.
   read(10,'(A)',iostat= ioerr) oneline
   if (ioerr.ne.0) then
      write(*,'(10x,"Error while reading the equilibrium density info from file: ",A)') trim(molfile)
      close(10)
      stop
   else
      read(oneline,*,iostat=ioerr) den0, TCdim
      if (ioerr.ne.0) then
         read(oneline,*,iostat=ioerr) den0
         if (ioerr.ne.0) then
            write(*,'(10x,"Error while reading the equilibrium density info from file: ",A)') trim(molfile)
            close(10)
            stop
         endif
         TCdim = 4
      endif
   endif
   !
   ! allocate memories
   allocate( EMass(ntypeMol), sigma(ntypeMol), epison(ntypeMol),           &
   &         bondcoef(2, nbndtypMol), anglcoef(2, nangtypMol),             &
   &         torscoef(TCdim, ndihtypMol), DistMin(ntypeMol), EName(ntypeMol+1) )
   allocate( bondsMol(3, nbondsMol), AnglsMol(4, nanglsMol),                &
   &         dihrsMol(5, ndihrsMol), attypMol(natomMol),                    &
   &         atposMol(3,natomMol), chargeMol(natomMol), molIDMol(natomMol) )                
   !
   ! The following lines are read by keywords, empty lines are skipped
   ! The availabe keywords are: Elemental Info, Bond Coeffs, Angle Coeffs, Dihedral Coeffs,
   ! Atoms, Bonds, Angles, Dihedrals
   read(10,'(A)',iostat=ioerr) oneline
   do while (ioerr.eq.0)
     if (oneline.eq.'Elemental Info') then
        ! Elementail info, each line carries the following info for each type of atom:
        ! id, name, mass (g/mol), sigma (Angstrom), epsilon (kcal/mol), dist-min
        do i = 1, ntypeMol
           read(10,'(A)', iostat=ioerr) oneline
           if (ioerr.ne.0) stop 'Error while reading elemental info from file!'
           read(oneline, *) id
           read(oneline, *) idum, EName(id), EMass(id), sigma(id), epison(id), DistMin(id)
        enddo
     elseif (oneline.eq.'Bond Coeffs') then
        ! Bond coefficient info, each line carries the coefficients for one type of bond
        ! id, K, r0
        do i = 1, nbndtypMol
           read(10,'(A)', iostat=ioerr) oneline
           if (ioerr.ne.0) stop 'Error while reading bond coeffs info from file!'
           read(oneline, *) id
           read(oneline, *) idum, bondcoef(:,id)
        enddo
     elseif (oneline.eq.'Angle Coeffs') then
        ! Angle coefficient info, each line carries the coefficients for one type of angle
        ! id, K, theta0
        do i = 1, nangtypMol
           read(10,'(A)', iostat=ioerr) oneline
           if (ioerr.ne.0) stop 'Error while reading angle coeffs info from file!'
           read(oneline, *) id
           read(oneline, *) idum, anglcoef(:,id)
        enddo
     elseif (oneline.eq.'Dihedral Coeffs') then
        ! Torsion coefficient info, each line carries the coefficients for one type of dihedral
        ! id, V1, V2, V3, V4
        do i = 1, ndihtypMol
           read(10,'(A)', iostat=ioerr) oneline
           if (ioerr.ne.0) stop 'Error while reading dihedral coeffs info from file!'
           read(oneline, *) id
           read(oneline, *) idum, torscoef(:,id)
        enddo
     elseif (oneline.eq.'Atoms') then
        ! Atom info, one line for one atom:
        ! id mol-id type charge x y z
        do i = 1, natomMol
           read(10,'(A)', iostat=ioerr) oneline
           if (ioerr.ne.0) stop 'Error while reading atoms info from file!'
           read(oneline, *) id
           read(oneline, *) idum, molIDMol(id), attypMol(id), chargeMol(id), atposMol(:,id)
        enddo
     elseif (oneline.eq.'Bonds') then
        ! Bonds info, one line for one bond
        ! id bond-type atom1 atom2
        do i = 1, nbondsMol
           read(10,'(A)', iostat=ioerr) oneline
           if (ioerr.ne.0) stop 'Error while reading bonds info from file!'
           read(oneline, *) id
           read(oneline, *) idum, bondsMol(:,id)
        enddo
     elseif (oneline.eq.'Angles') then
        ! Angles info, one line for one angle
        ! id angle-type atom1 atom2 atom3
        do i = 1, nanglsMol
           read(10,'(A)', iostat=ioerr) oneline
           if (ioerr.ne.0) stop 'Error while reading bonds info from file!'
           read(oneline, *) id
           read(oneline, *) idum, AnglsMol(:,id)
        enddo
     elseif (oneline.eq.'Dihedrals') then
        ! Dihedrals info, one line for one dihedral
        ! id dihedral-type atom1 atom2 atom3 atom4
        do i = 1, ndihrsMol
           read(10,'(A)', iostat=ioerr) oneline
           if (ioerr.ne.0) stop 'Error while reading bonds info from file!'
           read(oneline, *) id
           read(oneline, *) idum, DihrsMol(:,id)
        enddo

     endif

     read(10,'(A)',iostat=ioerr) oneline
   enddo
   close(10)
   !
   nMol      = maxval(molIDMol)
   molmass   = 0.D0
   do i = 1, natomMol
      molmass = molmass + EMass( attypMol(i) )
   enddo
   forall ( i=1:3 ) minbox(i) = maxval(atposMol(i,:)) - minval(atposMol(i,:))
   minlen = maxval( minbox )
   allocate( rotatedMol( 3, natomMol) )
   !
   subname = subnamep
return
end subroutine
